986 


JOURNAL  OF  MICROELECTROMECHANICAL  SYSTEMS,  VOL.  15,  NO.  4,  AUGUST  2006 


Modeling  the  Fluid  Dynamics  of  Electrowetting  on 

Dielectric  (EWOD) 

Shawn  W.  Walker  and  Benjamin  Shapiro,  Member,  IEEE 


Abstract — This  paper  discusses  the  modeling  and  simulation  of 
a  parallel-plate  Electrowetting  On  Dielectric  (EWOD)  device  that 
moves  fluid  droplets  through  surface  tension  effects.  We  model  the 
fluid  dynamics  by  using  Hele-Shaw  type  equations  with  a  focus  on 
including  the  relevant  boundary  phenomena.  Specifically,  we  show 
that  contact  angle  saturation  and  hysteresis  are  needed  to  predict 
the  correct  shape  and  time  scale  of  droplet  motion.  We  demonstrate 
this  by  comparing  our  simulation  to  experimental  data  for  a  split¬ 
ting  droplet.  Without  these  boundary  effects,  the  simulation  shows 
the  droplet  splitting  into  three  pieces  instead  of  two  and  the  motion 
is  over  15  times  faster  than  the  experiment.  We  then  show  how  in¬ 
cluding  the  saturation  characteristics  of  the  device,  and  a  simple 
model  of  contact  angle  hysteresis,  allows  the  simulation  to  better 
predict  the  splitting  experiment.  The  match  is  not  perfect  and  suf¬ 
fers  mainly  because  contact  line  pinning  is  not  included.  This  is 
followed  by  a  comparison  between  our  simulation,  whose  parame¬ 
ters  are  now  frozen,  and  a  new  experiment  involving  bulk  droplet 
motion.  Our  numerical  implementation  uses  the  level  set  method, 
is  fast,  and  is  being  used  to  design  algorithms  for  the  precise  con¬ 
trol  of  microdroplet  motion,  mixing,  and  splitting.  [1439] 

Index  Terms — Control,  electrowetting,  level  set  method,  mi¬ 
crofluidics,  modeling,  two-phase  flow. 


I.  Introduction 

WELL-DESIGNED  MEMS  devices  take  advantage  of  the 
large  surface  to  volume  ratios  found  at  the  microscale. 
In  particular,  microfluidic  devices  often  exploit  surface  tension 
forces  to  actuate  or  control  liquids  [1] — [3].  Electrowetting  refers 
to  using  electrical  fields  to  effectively  modify  surface  tension  ef¬ 
fects  [4]— [6]  (see  [7]  and  [8]  for  some  fascinating  experimental 
demonstrations). 

Applications  for  these  devices  range  from  microfluid  trans¬ 
port  [9],  mixing  [10],  dispensing  [11],  and  “lab-on-a-chip”  de¬ 
vices  that  automate  functions  like  sensing  and  testing  of  biolog¬ 
ical  samples  [12],  [13]  to  tunable  optical  fiber  devices  [14],  [15] 
and  reflective  displays  [16]. 

For  this  paper,  we  are  concerned  with  modeling  a  specific 
variant  of  electrowetting  called  Electrowetting-On-Dielectric 
(EWOD)  [17],  which  has  an  extra  insulating  layer  to  enhance 
its  operation.  See  [18]  for  an  initial  experimental  demonstration 
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and  [19]  for  an  analysis  of  the  advantage  of  using  a  dielectric 
insulating  layer  in  an  electrowetting  system. 

Similar  applications  exist  for  the  EWOD  device  as  well,  such 
as  mass  spectrometry  [20],  mixing  [21],  “lab-on-a-chip”  [22], 
micro-injection  [23],  and  particle  separation  and  concentration 
control  [24].  The  potential  uses  of  these  technologies  could  be 
for  controlled  mixing  of  chemicals  and  automated  DNA  testing. 

Ultimately,  these  applications  will  need  accurate  fluid  dynam¬ 
ical  control  in  order  to  execute  their  many  subtasks  (i.e.,  particle 
control,  precise  droplet  motion,  splitting,  optimal  mixing,  etc.). 
But  this  will  also  require  accurate  models  to  help  design  robust 
controllers  as  well  as  guide  device  optimization.  Furthermore, 
these  models  must  be  convenient  and  cheap  to  use  in  order  to 
fit  within  available  control  design  and  optimization  methodolo¬ 
gies. 

Other  modeling  efforts  of  EWOD  include  [17],  which  gives 
a  basic  model  of  how  the  device  parameters  affect  droplet  split¬ 
ting.  Equilibrium  models  for  the  shape  of  sessile  drops  on  a 
charging  dielectric  plate  are  given  in  [25]  and  [26].  In  particular, 
[25]  considers  a  conducting  liquid  on  top  of  an  insulating  layer 
and  the  effect  of  charge  trapping  at  high  voltage  on  contact  angle 
saturation.  In  [26],  it  is  shown  that  liquid  resistance  can  lead  to 
contact  angle  saturation  in  the  EWOD  devices.  An  alternative, 
lumped  parameter,  electromechanical  model  for  a  one-dimen¬ 
sional  (1-D)  liquid  column  actuated  by  electrowetting  is  given 
in  [27]  for  the  equilibrium  case  and  in  [28]  for  the  dynamics.  In 
addition,  a  dynamic  model  of  the  contact  angle  variation  for  a 
spreading  axisymmetric  drop  is  given  in  [29].  And  recently  in 
[30],  a  diffuse  interface  model  and  simulation  of  droplet  motion 
is  compared  to  experiments  on  a  scaled-up  version  of  the  elec¬ 
trowetting  device. 

In  this  paper,  we  present  a  distributed  parameter  model  of 
EWOD  fluid  dynamics  that  is  able  to  approximately  capture  the 
evolution  of  a  droplet’s  liquid-gas  interface,  using  the  level  set 
method  [3 1],  in  two  dimensions.  Our  model  includes  a  rough  ap¬ 
proximation  of  contact  angle  hysteresis,  which  is  different  than, 
though  analogous  to,  the  contact  line  friction  model  discussed  in 
[28]  and  [29].  Furthermore,  the  simulation  of  our  model  is  suf¬ 
ficiently  fast  and  low  dimensional  to  use  in  controller  design. 

This  paper  is  organized  as  follows.  Section  II  gives  an 
overview  of  the  EWOD  device’s  form  and  function.  Section  III 
develops  the  governing  fluid  equations  and  boundary  con¬ 
ditions,  along  with  our  model  of  contact  angle  hysteresis. 
Section  IV  describes  the  numerical  solution  scheme,  which 
uses  a  level  set  method  for  tracking  the  droplet  boundary.  And 
Section  V  presents  our  numerical  results  in  comparison  with 
two  experiments  that  exhibit  droplet  splitting  and  bulk  droplet 
motion. 
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Fig.  1.  Schematic  of  sample  EWOD  device  (courtesy  of  C.  J.  Kim  at  UCLA). 


I  |  Electrode  (Ground) 

Teflon 

I  I  Silicon  Dioxide 

HU  Electrode  (Grid) 

Fig.  2.  Cross-sectional  view  of  the  EWOD  device. 


II.  Description  of  the  Device 

A  schematic  of  an  EWOD  device  is  given  in  Fig.  1,  while 
Fig.  2  shows  a  cross-sectional  view.  The  device  consists  of  a 
sandwich  of  various  layers  listed  from  top  to  bottom  as:  top 
(transparent)  electrode,  hydrophobic  Teflon  coating,  droplets  of 
water  (here  only  one  droplet  is  shown),  another  Teflon  coating, 
a  layer  of  solid  dielectric  silicon  dioxide,  and  an  underlying  grid 
of  electrodes.  There  are  also  spacers  to  ensure  that  the  channel 
height  is  uniform. 

The  basic  principle  of  operation  is  that  the  liquid-gas  in¬ 
terface  of  the  droplet  can  be  locally  deformed  by  capacitively 
charging  the  silicon  dioxide  layer  underneath  it.  The  induced 
motion  of  the  droplet  is  due  to  competing  effects  of  energy 
storage  between  the  dielectric  layer  and  the  surface  energy  of  the 
liquid-gas  interface  [26].  Effectively,  each  electrode  can  change 
the  surface  tension  properties  immediately  above  it.  This  change 
can  be  used  to  move  droplets  from  electrode  to  electrode,  to  split 
droplets  (by  pulling  on  either  side  using  three  electrodes),  to  join 
droplets  by  making  them  collide,  and  to  mix  fluid  in  droplets  by 
making  the  droplets  execute  complex  paths. 

An  experimental  device  with  a  splitting  droplet  is  shown  in 
Fig.  3  (the  view  is  through  the  top  transparent  electrode).  The 
actuation  voltages  of  the  three  electrodes  from  left  to  right  are 
25,  0,  and  25  V. 

In  [26],  we  presented  a  model  for  the  equilibrium  shape  of 
droplets  under  applied  electric  fields.  In  this  paper,  we  further 
consider  the  nonequilibrium  fluid  dynamics.  Specifically,  we 
focus  on  modeling  and  simulating  motion,  splitting,  and  joining 
of  the  liquid  droplets. 

III.  EWOD  Modeling 

This  section  describes  the  EWOD  modeling  approach.  In  par¬ 
ticular,  our  main  assumptions,  derivation  of  the  fluid  equations, 


Fig.  3.  Top  view  of  experimental  EWOD  device  (courtesy  of  C.  J.  Kim  at 
UCLA). 


TABLE  I 

Physical  Parameters 

Parameter  Symbol  Definition 

Surface  Tension  <r;s  =  0.07199  J/m' 

Dynamic  Viscosity  p  =  0.89  g/m  ■  s 
Density  p  =  996.93  Kg/m3 

Channel  Height  H  =  70  pm 

Electrode  Length  Lsiec  =  1.4  mm 

Length  Scale  L«3x  Leu-c 

Velocity  Scale  U0  (see  Sec.  V-A) 

Time  Scale  to  =  L/Uo 

Pressure  Scale  Po  =  <Jig/L 

Reynolds  Number  Re  =  pUoH/p 

Capillary  Number  Ca  =  pUo/aig 


Fig.  4.  EWOD  device  geometry.  The  coordinate  axes  are  defined  such  that  the 
top  and  bottom  plates  of  the  device  lie  in  planes  parallel  to  the  g-y  plane.  The 
physical  parameters  of  the  device  are  listed  in  Table  I. 


proper  boundary  conditions,  voltage  actuation,  contact  angle 
saturation,  and  hysteresis  effects  are  discussed  in  detail.  A  list  of 
the  physical  parameters  for  the  geometry  of  the  EWOD  device, 
as  well  as  the  fluid  parameters  for  distilled  water  at  standard 
temperature  and  pressure  (assumed  in  our  model),  is  given  in 
Table  I. 

A.  Governing  Equations  of  the  Liquid  Flow 

In  the  following  sections,  the  main  assumptions  and  gov¬ 
erning  equations  for  the  flow  of  liquid  inside  an  EWOD  device 
are  described  (see  Fig.  4).  In  particular,  we  obtain  a  model  sim¬ 
ilar  to  Hele-Shaw  type  flow  with  pressure  boundary  conditions 
at  the  liquid-gas  interface  proportional  to  its  mean  curvature. 

1)  Navier-Stokes  Equations:  We  start  by  considering  the 
Knudsen  number  of  the  EWOD  device,  which  provides  a  mea¬ 
sure  of  how  accurate  the  continuum  hypothesis  is  for  a  fluid 
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system  [32].  For  our  case,  we  can  assume  the  flow  physics  to 
be  a  continuum  because  the  Knudsen  number  is 


Kn  = 


H 


6.111  x  nr8  m 

70  x  10-6  m 


=  8.73  x  1(T4 


where  Aa;r  is  the  mean  free  path  of  air  molecules  at  standard 
temperature  and  pressure,  and  H  is  the  channel  height  of  the 
device.  Clearly,  this  is  within  the  continuum  regime  which  is 
defined  to  be  Kn  <  10-2.  The  Knudsen  number  associated  with 
the  liquid  flow  is  even  smaller  because  the  mean  free  path  of 
water  is  much  lower  than  that  of  air. 

Since  the  flow  is  a  continuum,  the  dimensional  Navier-Stokes 
equations  are  applicable.  Because  we  are  modeling  the  flow  of 
water,  incompressibility  and  Newtonian  fluid  assumptions  may 
be  used  [33].  This  gives 


V  •  V  =  0  (1) 

p(Vt  +  V  ■  VV)  =  -VP  +  p,V2V  (2) 

in  the  bulk  liquid,  where  V  =  (u.  v,  w )  is  the  three  dimensional 
velocity,  P  is  the  pressure,  subscript  t  denotes  the  partial  deriva¬ 
tive  with  respect  to  time,  and  p  and  //  are  the  density  and  dy¬ 
namic  viscosity,  respectively.  Equations  (1)  and  (2)  represent 
conservation  of  mass  and  momentum,  respectively,  with  gravity 
ignored  because  the  potential  energy  change  in  the  z  direction 
is  negligible  when  the  channel  height,  H,  is  small. 

Next,  we  have  the  boundary  conditions  for  a  liquid  droplet 
between  two  parallel  plates.  On  the  top  and  bottom  plates,  we 
have  the  usual  no-slip  condition  for  velocity  (i.e.,  all  velocity 
components  are  zero).  Since  the  air  surrounding  the  droplet  is 
not  being  forced,  it  does  not  significantly  affect  any  droplet  mo¬ 
tion.  Therefore,  by  ignoring  the  airflow,  we  have  the  following 
conditions  for  the  free  surface  of  an  incompressible,  Newtonian 
liquid  (i.e.,  the  liquid-gas  interface)  [34] 

n  ■  T n  =  — cr\gK  (3) 

t  ■  Tn  =  0  (4) 

where  erig  denotes  surface  tension,  k  is  the  mean  curvature  of 
the  interface  [35],  T  is  the  stress  tensor,  h  is  the  unit  normal 
vector  to  the  interface,  and  t  is  any  tangent  vector  to  the  in¬ 
terface.  Physically,  (3)  states  that  the  normal  stress  across  the 
liquid-gas  interface  is  balanced  by  surface  tension,  whereas  (4) 
says  the  tangential  stress  vanishes  because  the  airflow  is  negli¬ 
gible, 

2)  Hele-Shaw  Type  Flow:  Since  we  have  pressure-driven 
flow  in  a  slot  with  channel  height  much  smaller  than  the  diam¬ 
eter  of  the  droplet  [33],  the  Reynolds  number  is  small  and  we 
assume  the  flow  can  be  modeled  in  two  dimensions.  By  making 
the  additional  assumption  that  the  x  and  y  fluid  velocity  com¬ 
ponents  u  and  v  have  a  quadratic  profile  in  the  z  direction  (i.e., 
local  Poiseuille  flow;  see  Fig.  5),  (1)  and  (2)  can  be  nondimen- 
sionalized  and  reduced  to  a  form  similar  to  Hele— Shaw  flow 
[34]; 

V2P  =  0  (5) 

(»“‘  +  12(s)2"=-^  (6) 


Fig.  5.  Velocity  profile:  the  fluid  velocity  field  is  assumed  to  have  a  quadratic 
profile  in  the  0  direction. 


Fig.  6.  Overhead  view  of  a  (2-D)  droplet  with  side  view  zoom-in  of  the  inter¬ 
face.  The  liquid-gas  interface  is  assumed  to  have  a  circular  cross-section,  which 
gives  an  estimate  of  the  z  curvature,  nz,  in  dimensional  form.  The  x-y  curva¬ 
ture,  nxy,  is  just  the  curvature  of  the  boundary  of  the  2-D  droplet. 


(H-Ki)’— (7) 

where  the  term  on  the  far  left  of  (6)  and  (7)  is  the  extra  term  be¬ 
yond  the  usual  Hele-Shaw  equations.  This  time  derivative  term 
is  included  because  it  may  have  a  large  magnitude  due  to  rapidly 
varying  pressure  boundary  conditions  if  high-frequency  voltage 
actuation  is  used  to  modulate  the  droplet’s  contact  angles. 

The  boundary  conditions  for  (5)  are  then  given  by  the 
Young-Laplace  relation  [34],  which  says  (in  nondimensional 
form)  that  the  pressure  on  the  liquid-gas  interface  is  equal 
to  the  mean  curvature  of  the  interface.  Because  the  channel 
spacing  is  so  small,  this  is  accurately  approximated  by 

P  =  Kxy  +  yyKzi  at  the  liquid/gas  interface  (8) 
12 

where  Kxy  is  the  nondimensional  curvature  of  the  droplet  in  the 
x-y  plane,  kz  is  the  nondimensional  curvature  of  a  cross-sec¬ 
tion  of  the  droplet  along  the  z  axis  (see  Fig.  6),  and  L  is  the 
x-y  length  scale  of  the  device.  Given  that  (5 )  has  been  posed  in 
two  dimensions,  (8)  is  evaluated  at  each  point  of  the  two-dimen¬ 
sional  (2-D)  droplet  boundary  and  is  discussed  in  Section  III-B 1 . 


B.  Physics  of  the  Droplet  Boundary 

Above,  we  described  the  governing  equations  of  liquid 
droplet  motion.  We  now  discuss  the  geometry  and  different 
physical  phenomena  happening  at  the  liquid-gas  interface,  such 
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Fig.  7.  Curvature  note:  k xy  and  k  are  both  positive  for  the  bulging  droplet  on 

top.  For  the  inward  bending  droplet,  .  is  negative.  where  the  coefficients  are  given  by 


as  voltage  actuation,  contact  angle  saturation,  and  hysteresis. 
We  show  how  these  effects  are  modeled  and  how  they  affect 
the  computation  of  the  boundary  conditions. 

1 )  Interface  Curvature:  The  interface  mean  curvature  is  ap¬ 
proximated  using  the  individual  curvatures  nxy  and  nz  in  (8). 
We  compute  the  2  curvature  by  assuming  the  interface  has  a 
circular  cross-section  (see  Fig.  7).  The  x-y  curvature  compu¬ 
tation  requires  a  representation  of  the  shape  of  the  two  dimen¬ 
sional  droplet  boundary.  This  is  accomplished  by  using  a  level 
set  method  to  implicitly  capture  the  interface  and  is  described 
in  more  detail  in  Section  IV-C2. 

To  use  the  circular  approximation  for  computing  the  2  curva¬ 
ture,  we  must  know  the  slope  of  the  liquid-gas  interface  cross- 
section  at  the  floor  and  ceiling  of  the  EWOD  device.  This  is 
given  by  the  top  and  bottom  contact  angles,  0t  and  @b  respec¬ 
tively  (see  Fig.  6).  After  some  basic  geometry,  the  dimensional 
z  curvature  is  given  by 


Kz 


1 

H 


[cos(flt)  +  cos(0&)] 


which  gives  the  nondimensional  curvature  as 


Kz  =  -[cos(0t)  +  cos(06)]  (9) 

Through  the  use  of  EWOD,  contact  angles  can  be  voltage 
controlled,  thereby  allowing  the  droplet  pressure  field  to  be 
actuated  through  the  boundary  condition  (8).  The  details  of 
voltage  actuation  and  contact  angles  are  described  in  the  next 
Sections  II1-B2  and  III-B3. 

2)  EWOD  Charging  Time:  We  analyze  the  electrical 
charging  time  of  the  EWOD  device  to  determine  if  it  must  be 
included  in  our  model.  Consider  the  circuit  diagram  shown 
in  Fig.  8.  Using  transfer  function  theory  [36],  which  reduces 
the  ordinary  differential  equation  associated  with  the  electric 
circuit  to  an  algebraic  problem,  we  can  obtain  an  estimate  for 
the  time  constant  of  the  device.  The  transfer  function  for  this 
circuit  is  given  by 


T(s) 


ais  +  CK2 
0:3s2  +  0:4s  +  0:5 


(10) 


Q-i  =  (Ct  +  Co)RoRt 

Ot-2  =  Ro  +  Rt 

<23  =  RlRoCoRtCt 

0:4  =  RoRt(Ct  +  2Co)  +  Rl(RoCo  +  RtCt) 
CK5  =  Ro  +  Rl  +  2  Rt- 


Using  the  parameters  for  water,  silicon  dioxide.  Teflon,  and 
the  dimensions  of  the  EWOD  device  [17],  (10)  becomes 

1927.5(s  +  3.957  x  10_6) 

“  0  +  3.132  x  10-6)(s  + 2439.7)' 

After  approximately  canceling  the  two  near-identical  terms  in 
the  numerator  and  denominator,  we  are  left  with  a  transfer  func¬ 
tion  describing  a  first-order  differential  equation.  The  defining 
parameter  of  any  first  order  system  is  the  time  constant,  which 
in  this  case  is  0.41  ms.  Using  this,  the  electrical  charging  time 
is  estimated  as  four  times  the  time  constant,  or  1.64  ms. 

For  the  splitting  droplet  experiment  in  Section  V-A,  the  time 
to  split  is  0.167  s.  Since  the  majority  of  the  voltage  drop  occurs 
across  the  bottom  SiO_>  and  Teflon  layer,  and  the  charging  time 
is  over  100  times  faster  than  the  bulk  fluid  motion  we  are  inter¬ 
ested  in  (i.e.,  droplet  splitting),  we  assume  the  output  voltage 
is  instantaneously  equal  to  the  input  voltage.  Therefore,  given 
that  there  is  a  direct  relation  between  contact  angle  and  applied 
voltage  (see  Section  III-B3),  EWOD  is  capable  of  changing  the 
contact  angle  very  quickly. 

3 )  Contact  Angles  and  Saturation:  There  is  a  considerable 
amount  of  literature  on  contact  angles  and  wetting  phenomena; 
see  the  following  references  for  a  sampling  [37]— [40] .  For  a 
detailed  physical  description  of  electrowetting,  see  [41]— [44], 
[25],  In  this  section,  we  are  concerned  with  how  the  contact 
angle  varies  with  respect  to  the  applied  voltage. 

According  to  [4],  [45],  and  [29],  for  a  sessile  drop  on  a  single 
dielectric  plate,  the  Young-Lippmann  equation  predicts  a  near 
parabolic  curve  relating  contact  angle  to  the  capacitive  voltage 
across  the  plate  (see  Fig.  9).  However,  if  Young — Lippmann  is 
used  to  simulate  droplet  splitting,  it  predicts  an  incorrect  shape 
for  the  motion  of  the  droplet.  This  is  because  electrowetting, 
in  reality,  deviates  from  the  Young-Lippmann  theory  at  high 
voltages  and  reaches  a  saturation  limit  (also  shown  in  Fig.  9). 
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Contact  Angle  Vs.  Voltage 


Fig.  9.  Contact  angle  versus  voltage  curves.  Theoretical  and  experimental  data 
for  contact  angle  variations  using  electrowetting  on  dielectric  are  presented.  The 
dotted  line  denoting  the  Young-Lippmann  (Y-L)  curve  is  theoretical  [4],  [45], 
[29].  The  single  plate  saturation  curve  has  six  experimental  data  points  (given 
in  [17])  with  a  piecewise  interpolating  polynomial  (see  dashed  line  and  “o” 
data  points).  The  two  plate  saturation  curve  has  two  experimental  data  points 
(also  from  [17])  with  an  interpolating  curve  derived  from  the  single  plate  case 
in  Section  III  using  a  linear  map  (see  solid  line  and  “o”  data  points).  This  curve 
is  used  in  our  simulation  to  predict  the  correct  droplet  motion  and  splitting  time 
(see  Section  V). 


In  Sections  V-Al  and  V-A2,  we  present  simulations  using  the 
Young-Lippmann  equation  and  saturation,  respectively,  to  illus¬ 
trate  the  importance  of  modeling  the  latter.  For  more  informa¬ 
tion  on  the  causes  of  contact  angle  saturation  of  sessile  droplets, 
see  [26], 

The  available  literature  only  discusses  a  sessile  droplet  on 
a  single  plate.  For  this  paper,  we  need  data  on  contact  angle 
variations  of  a  droplet  sandwiched  between  two  plates.  In  [17], 
experimental  contact  angle  data  for  the  EWOD  device  is  given  at 
an  applied  voltage  of  0  and  25  V.  The  top  contact  angle  remains 
approximately  the  same  at  117°  regardless  of  voltage  actuation. 
This  is  because  most  of  the  dielectric  energy  is  stored  in  the 
bottom  plate  due  to  the  presence  of  the  SiCD  layer.  Therefore, 
we  assume  the  nominal  contact  angle  on  the  top  plate  is  fixed  at 
117°.  The  bottom  contact  angle  varies  between  117°  and  90°  at 
0  and  25  V,  respectively. 

In  order  to  model  contact  angle  variations  on  the  bottom  plate 
for  any  voltage,  we  must  combine  the  two  experimental  data 
points  for  the  parallel  plate  case  with  the  six  data  points  for  the 
single  plate  case  (shown  in  Fig.  9).  Since  there  is  an  overlap 
between  two  of  the  data  points  for  both  cases,  we  define  a  linear 
mapping  that  represents  this  correspondence 


1170  _  qno 

02  =— - ^V(01  -80.4°)  +90°  (11) 

119°  -  80.4° v  ' 

where  6i  is  the  contact  angle  for  a  single  plate  and  6o  is  the 
contact  angle  for  two  plates.  This  equation  maps  1 19.0°  to  117° 
and  80.4°  to  90°.  By  taking  the  six  data  points,  and  their  piece- 
wise  interpolating  polynomial,  for  saturation  on  a  single  plate 
as  input  to  this  linear  map,  we  obtain  the  saturation  curve  for 
two  plates  as  output  (see  Fig.  9).  Due  to  the  scarcity  of  data  on 


I 

9 

-\ 

Fig.  10.  Illustration  of  contact  line  pinning  and  hysteresis.  The  contact  line 
of  the  droplet  is  where  the  liquid-gas  interface  meets  the  solid  surface.  Line 
pinning  simply  means  the  contact  line  (and  the  droplet)  is  stuck  to  the  surface. 
A  direct  result  of  this  is  contact  angle  hysteresis,  which  refers  to  the  situation 
where  the  receding  angle  9r  is  less  than  the  nominal  (equilibrium)  angle  9$, 
while  the  advancing  angle  9  J 4  is  greater  than  90 .  In  the  diagram  above,  90  is  the 
contact  angle  of  the  droplet  on  a  horizontal  surface,  whereas  9r  and  9  a  are  the 
contact  angles  when  the  surface  is  tilted.  The  droplet  can  slide  by  using  a  large 
enough  tilt  angle,  but  the  motion  will  be  limited  by  the  static  frictional  effect  of 
line  pinning  and  contact  angle  hysteresis  will  still  be  present.  A  similar  situation 
happens  in  EWOD,  where  hysteresis  also  acts  as  a  retarding  effect  by  deforming 
the  liquid-gas  interface  shape  in  an  unfavorable  way  (see  Section  III-B4). 


contact  angle  variation  for  the  parallel  plate  EWOD  device,  we 
assume  the  two  plate  saturation  curve  in  Fig.  9  is  true  for  our 
model. 

4)  Hysteresis:  Contact  angle  hysteresis  is  the  last  piece  of 
physics  we  need  to  complete  our  model  of  droplet  motion  using 
EWOD  forces.  Hysteresis  refers  to  the  difference  in  contact  an¬ 
gles  between  the  advancing  and  receding  ends  of  sessile  drops. 
It  is  a  direct  consequence  of  contact  line  pinning,  which  acts  as 
a  force  that  resists  any  sliding  motion,  and  can  be  seen  when 
water  droplets  stick  to  the  side  of  a  solid  surface  (see  Fig.  10). 
For  more  information  on  contact  angle  hysteresis  and  line  pin¬ 
ning,  see  [37],  [39],  and  [40], 

From  Fig.  10,  for  a  sessile  drop  on  a  single  plate,  it  can  be 
seen  that  the  advancing  and  receding  contact  angles  are  greater 
and  smaller,  respectively,  than  the  nominal  contact  angle.  This 
is  also  true  for  a  droplet  inside  the  EWOD  device  (shown  in 
Fig.  11).  Ideally,  if  there  were  no  hysteresis,  the  nominal  contact 
angle  at  the  interface  of  the  droplet  should  be  determined  by  the 
two  plate  saturation  curve  in  Fig.  9  and  the  applied  voltage  at  the 
interface.  But  in  the  presence  of  hysteresis,  the  contact  angles 
deflect  from  their  nominal  values  which  affects  the  pressure  on 
the  boundary  by  the  Young-Laplace  relation  (8). 

To  see  how  it  is  affected,  consider  a  circular  droplet  in  mo¬ 
tion  due  to  voltage  actuation  (see  Fig.  11).  Let  PA  and  Pr  de¬ 
note  the  pressures  at  the  advancing  and  receding  ends  of  the 
droplet,  respectively,  when  no  hysteresis  is  present.  And  let  Pa 
and  Pj i  denote  the  same  pressures  with  hysteresis.  It  is  clear 
from  Fig.  1 1  that  the  z  curvatures  at  the  receding  and  advancing 
ends  of  the  droplet  are  larger  and  smaller,  respectively,  for  no 
hysteresis  than  with  hysteresis.  From  (8),  it  can  be  seen  that  the 
curvature  change  implies  that  Pr  >  Pr  and  Pa  <  PA. 

This  change  in  boundary  pressure  weakens  the  pressure  gra¬ 
dient  throughout  the  droplet  from  what  it  would  be  without  hys¬ 
teresis  because  its  magnitude  is  proportional  to  the  pressure  dif¬ 
ference 


|VP|  oc  \Pr  —  PA\ 
|VP|  oc  \PR  —  PA\ 
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Pr  >  Pr  With  Hysteresis  Pa  >  Pa 


OFF  |  ON 


Fig.  11.  Effects  of  contact  angle  hysteresis  in  the  EWOD  device.  A  droplet 
is  shown  moving  from  left  to  right  due  to  voltage  actuation  (OFF/ON).  When 
hysteresis  is  present,  the  contact  angles  differ  from  their  nominal  (nonhysteresis) 
values.  The  effects  on  the  liquid-gas  interface  pressure  are  also  shown.  The 
presence  of  hysteresis  causes  the  pressure  gradient  throughout  the  droplet  to  be 
weakened  from  the  nominal  case  (Pr  —  Pa  >  Pr  —  Pa)- 


where  VP  and  VP  are  the  pressure  gradients  with  and  without 
hysteresis,  respectively.  Using  the  fact  that  Pr  >  Pr  >  P4  > 
P4,  we  obtain  the  inequality 


|VP|  <  |VP|. 


where  the  pressure  terms  are  given  by  (8)  and  (9).  Using  these, 
ivhys  can  be  estimated  by 


U'i  mm 
1-Men 
^Mtys 


[cos(6poy  Ahys)  T  oosftpov  Ahys) 

-  COs(0t,25V  +  Ahys)  -  COs(6p25V  +  Ahys)] 

[cos(0tfov)  +  cos(8b,ov) 

-  cos(0t,25v)  -  COs(0fc,25v)] 


Gn 


1-Men 


(12) 


where  Ahys  is  the  extra  amount  of  contact  angle  deflection  from 
the  nominal  angle  due  to  hysteresis,  and  8t  and  8b  are  the  top 
and  bottom  contact  angles,  respectively.  The  voltage  subscripts 
specify  the  actuation  strength  of  the  contact  angles,  with  the  top 
angle  fixed  at  117°  regardless  of  voltage  and  the  bottom  angle 
obeys  the  two  plate  saturation  curve  in  Fig.  9.  Due  to  the  lack  of 
data  on  hysteresis  of  EWOD  driven  droplets,  we  assume  the  con¬ 
tact  angle  deflection  to  be  the  same  on  the  top  and  bottom  of  the 
advancing  and  receding  ends  of  the  droplet.  In  Section  V.A.3, 
we  use  (12)  to  estimate  the  hysteresis  angle  deflection  that  cor¬ 
responds  to  the  appropriate  constant,  K^ys,  that  ensures  the  sim¬ 
ulated  splitting  time  matches  the  experiment. 


Hence,  the  driving  force  of  the  droplet  motion  is  decreased  when 
hysteresis  is  present.  This  is  why  our  simulation  (with  just  satu¬ 
ration  modeled)  predicts  a  split  time  over  nine  times  faster  than 
the  experiment  shows  (see  Section  V.A). 

From  the  discussion  above,  an  obvious  way  to  model  hys¬ 
teresis  is  to  modify  the  contact  angle  of  the  interface  based  on 
which  way  it  is  moving.  However,  from  our  own  numerical  ex¬ 
periments,  this  method  is  not  very  robust.  Therefore,  we  opted 
for  a  simpler  model  by  assuming  that 

VP  =  KhysVP 

where  K\ys  is  a  constant  smaller  than  one.  In  other  words,  we 
scale  down  the  pressure  gradient  in  (6)  and  (7)  to  account  for 
hysteresis.  This  is  analogous  to  the  contact  line  friction  model 
in  [29]  and  [28],  which  also  acts  as  a  retarding  effect  on  liquid 
motion. 

Scaling  the  pressure  gradient  is  simple  and  straightforward, 
yet  is  capable  of  approximately  capturing  the  droplet  motion 
and  time  scale  observed  in  the  experiments  (see  Section  V.A.3). 
However,  we  do  stress  that  it  is  not  exact.  It  does  not  capture  the 
effect  of  line  pinning,  which  is  observable  in  EWOD  as  demon¬ 
strated  in  [  17]  by  the  fact  that  droplets  do  not  move  unless  suf¬ 
ficient  voltage  actuation  is  used.  Line  pinning  is  not  completely 
understood  and  implementing  a  highly  accurate  model  would 
be  computationally  expensive.  Hence,  we  opted  for  a  model 
that  is  computationally  quick  but  still  captures  the  lossy  effect 
of  droplet  motion  that  is  induced  by  line  pinning;  namely  hys¬ 
teresis. 

We  now  estimate  the  hysteresis  constant  in  terms  of  contact 
angles.  From  the  relations  given  above,  we  have 


C.  Final  Equation  Summary 

We  now  write  the  final  model  equations  describing  the  fluid 
flow  of  a  liquid  droplet  inside  an  EWOD  device.  The  equations 
for  the  pressure  field  are 

V2P  =  0,  in  n  (13) 

P  =  Kxy  +  on  OQ  (14) 

11 

where  V2  is  the  Laplacian  operator,  Q  denotes  the  domain  of 
the  liquid  droplet  in  two  dimensions,  dfl  is  its  boundary,  P  is 
the  pressure,  L  is  a  chosen  length  scale,  H  is  the  channel  height, 
and  nxy  is  the  curvature  in  the  x-y  plane.  The  2  curvature,  kz, 
is  given  by 

=  -[cos  (9t)  +  cos  (6b)]  (15) 

where  8t  and  8b  are  the  contact  angles  on  the  top  and  bottom  of 
the  EWOD  device,  respectively.  The  top  angle  is  assumed  to  be 
117  degrees  regardless  of  the  applied  voltage.  The  variations  of 
the  bottom  angle  are  given  by  the  two  plate  saturation  curve  in 
Fig.  9. 

The  equations  for  the  velocity  field  are 

TUt  +  I3u  =  in^  (16) 

TVt  +  f3v  =  in^  (17) 

where  u  and  v  are  the  velocity  components  in  the  x  and  y  di¬ 
rections,  /lhys  is  tlle  hysteresis  constant,  and  Ca  is  the  capillary 
number.  The  constants  r  and  fl  are  given  by 
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(18) 
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where  Uo  is  a  chosen  velocity  scale,  p  is  the  fluid  density,  and  p 
is  the  dynamic  viscosity. 

Because  we  have  two-phase  flow,  we  need  the  following 
equation  to  describe  the  time-varying  nature  of  the  droplet 
domain 


dh  =  n  ■  (u,v)\an  (19) 

where  n  is  the  unit  normal  vector  of  the  boundary.  Basically,  the 
droplet  boundary  moves  with  the  normal  velocity  of  the  fluid.  In 
the  next  section,  we  discuss  the  numerical  simulation  of  these 
equations. 


IV.  Numerical  Implementation 

The  most  crucial  part  of  simulating  the  fluid  equations  in 
Section  III-C  is  in  handling  the  time-varying  two-phase  droplet 
boundary  modeled  by  (19).  Various  methods  for  simulating 
two-phase  flow  are  available  (see  [46J  for  a  survey).  For  this 
paper,  the  method  needs  to  be  capable  of  simulating  splitting 
and  joining  of  droplets  without  excessive  computation.  There¬ 
fore,  we  chose  the  level  set  method  (see  [31],  [47],  [48]),  which 
implicitly  tracks  the  droplet  boundary  as  the  zero  level  set  of  a 
scalar  function  defined  over  the  x-y  plane  (see  Fig.  14).  This 
scalar  function  deforms  and  changes  shape  based  on  the  fluid 
velocity  field,  which  causes  the  zero  level  set  to  also  change. 
Hence,  the  motion  of  the  droplet  boundary  is  captured  through 
the  evolution  of  the  scalar  function. 

This  level  set  function  is  used  to  define  the  domain  of  the 
droplet  at  each  instant  of  time,  allowing  the  pressure  and  ve¬ 
locity  fields  to  be  computed  from  finite  difference  approxima¬ 
tions  to  (13),  (14),  (15),  (16),  (17).  We  combine  these  methods  in 
a  third  order  Runge-Kutta  time-stepping  algorithm  which  spec¬ 
ifies  an  order  to  the  computation  of  the  pressure  field,  velocity 
field,  and  level  set  update  (see  Fig.  12).  The  following  sections 
give  the  details  of  our  algorithm  which  is  based  on  the  methods 
used  in  [49]  and  [48], 

All  simulations  were  performed  with  MATLAB  on  a  Pen¬ 
tium  4,  3.6  GHz  with  2  GB  of  RAM  running  Windows  XP.  The 
computing  time  of  each  simulation  varied  between  three  and  six 
minutes  for  a  108  x  108  mesh,  which  shows  the  speed  of  our 
method. 

A.  Discretization 

The  computational  domain  is  defined  to  be  the  unit  square, 
and  is  discretized  by  a  regular  cartesian  grid  (see  Fig.  13).  For 
the  simulations  given  in  Section  V,  we  used  a  108  X  108  mesh. 
On  this  grid,  the  level  set  function,  </>,  and  the  fluid  variables 
u,  v,  and  P  are  sampled.  A  small  buffer  region,  three  grid  nodes 
thick,  is  defined  at  the  sides  of  the  computational  domain.  No 
droplet  motion  is  allowed  inside  the  buffer  region  because  of 
potential  problems  with  computing  second  order  spatial  deriva¬ 
tives  there. 

B.  Initialization 

The  level  set  function,  <p,  is  initialized  to  a  signed  distance 
function  with  the  zero  level  contour  corresponding  to  the  initial 
interface  shape  (see  Fig.  14).  By  distance  function,  we  mean 


Correct  Level  Set 
Choose  Time  Step 


Fig.  12.  Algorithm  flowchart. 


0  1 


•  -  Interior  Nodes 
o  -  Boundary  Nodes 


Fig.  13.  Computational  domain  layout  (liquid  region  corresponds  to  interior 
nodes). 


that  the  value  of  (f>  at  a  grid  point  in  the  computational  domain 
corresponds  to  the  shortest  distance  that  the  grid  point  is  from 
the  interface.  Signed  distance  means  that  d>  is  positive  inside  the 
droplet  and  negative  outside.  Next,  the  velocity  field,  (u,v),  is 
set  to  zero.  And  finally,  we  choose  a  small  initial  time  step  before 
entering  the  main  update  routine  discussed  in  Section  IV-C. 

C.  Main  Update  Routine 

At  each  time  step  of  our  simulation,  the  fluid  variables  and 
level  set  function  are  updated  by  computing  a  convex  combina¬ 
tion  of  three  forward  Euler  steps.  This  method  is  a  third  order 
Runge-Kutta  method,  and  is  detailed  in  [31]  and  [50]. 

In  each  Euler  step,  the  level  set  is  updated  first,  followed  by 
the  pressure,  and  then  velocity.  The  updated  level  set  is  used  in 
computing  the  pressure  field  for  the  new  time  step,  which  is  then 
used  to  update  the  velocity  field  (see  Fig.  12).  In  the  following 
sections,  we  give  the  details  of  each  of  these  subroutines. 
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Fig.  14.  Example  of  a  level  set  function.  The  zero  level  contour  is  highlighted 
with  a  thick  black  line  and  shows  a  droplet  about  to  split. 


1)  Update  Level  Set:  The  level  set  function  is  updated 
through  a  standard  convection  equation: 


(ft  +  (u,v)  ■  V0  =  0  (20) 

This  equation  represents  conservation  of  the  quantity,  <f>,  while 
being  transported  by  the  velocity  field,  (u,v).  Hence,  the  zero 
level  set  “quantity,”  f  =  0,  is  transported  by  the  local  fluid  ve¬ 
locity  around  the  droplet  boundary.  The  numerical  implemen¬ 
tation  uses  an  upwinded  third  order  Hamilton-Jacobi  weighted 
essentially  nonoscillatory  (WENO)  method  for  discretizing  the 
velocity-gradient  term,  which  uses  u,  v,  and  from  the  pre¬ 
vious  time  step.  This  method  is  robust  and  is  described  in  detail 
in  [31], 

In  this  paper,  to  reduce  simulation  time,  only  the  grid  nodes  in 
a  small  band  surrounding  the  zero  level  set  are  actually  updated. 
This  does  not  reduce  accuracy  since  the  whole  level  set  must  be 
reset  to  a  distance  function  periodically  (see  Section  IV-D1). 

2)  Solve  for  Pressure:  Next,  the  domain,  f l,  of  the  droplet 
is  defined  to  be  the  regions  in  the  x-y  plane  where  the  level  set 
function,  <f>,  is  positive  (see  Fig.  14).  The  computational  domain 
of  a  hypothetical  droplet  is  depicted  in  Fig.  13.  Each  of  the  grid 
nodes  is  located  on  an  electrode  with  a  known  applied  voltage. 
The  local  x-y  curvature  of  the  boundary.  Oil,  is  then  given  by 
131]: 


4*x 4>H<I  2 4>x4>y4>xy  3”  4>y4)xx 


(« 4>l  +  4>l) 


a/2 


(21) 


b= 0 


where  the  level  set  derivative  terms  are  approximated  using  cen¬ 
tral  differences.  Because  of  fundamental  problems  with  differ¬ 
entiating  numerical  data,  the  level  set  function,  <f>,  must  be  fil¬ 
tered  prior  to  computing  the  derivative  terms  [50].  In  addition, 
Kxy  must  be  post-filtered  to  ensure  smooth  curvature  data.  This 


is  mainly  due  to  the  explicit  nature  of  the  curvature  calculation 
used  here. 

Then,  we  get  the  bottom  contact  angle,  6b,  at  each  boundary 
node  using  the  known  voltage  there  and  the  two  plate  satura¬ 
tion  curve  in  Fig.  9.  Voltage  transitions  near  the  edge  between 
two  electrodes  are  smoothed  out  using  linear  interpolation,  in 
a  narrow  region,  to  prevent  large  localized  velocities  caused  by 
discontinuous  boundary  conditions.  Finally,  the  boundary  pres¬ 
sure  values  are  computed  using  (14),  (15),  and  (21)  evaluated 
on  the  boundary  nodes. 

The  pressure  values  at  the  interior  nodes  are  computed  by 
solving  (13),  which  implicitly  contains  the  conservation  of  mass 
( 1 ).  The  numerical  solution  is  obtained  by  using  a  simple  Gauss- 
Seidel  iterative  solver  with  a  relative  error  tolerance  of  10-8 
[51],  Other,  more  advanced  methods  for  solving  a  matrix  system 
of  equations  exist,  but  would  require  the  matrix  structure  to  be 
recreated  at  every  time  step  since  the  domain  of  the  droplet  is 
always  changing.  In  addition,  the  Gauss-Seidel  solver  is  imple¬ 
mented  in  C,  for  speed,  and  called  from  MATLAB.  Therefore, 
we  saw  no  significant  advantage  with  using  a  different  method. 

Once  the  pressure  values  are  known,  the  pressure  gradient, 
( Px,Py ),  at  every  interior  node  is  computed  using  a  central  dif¬ 
ference  formula  [51].  These  values  are  then  used  in  the  velocity 
update  routine. 

3)  Update  Velocity  Field:  The  fluid  velocity  components, 
(u,  v),  obey  two  first  order  time  differential  equations  given  by 
(16)  and  (17).  The  pressure  gradient  provides  a  forcing  term  in 
the  equations,  which  causes  a  velocity  field  to  develop.  We  com¬ 
pute  the  velocities  on  our  computational  domain  by  discretizing 
(16)  and  (17)  in  space  while  keeping  time  continuous.  This  ap¬ 
proach  is  commonly  known  as  a  semi-discrete  method  [52]  (or 
method-of-lines)  and  allows  for  the  use  of  an  analytic  solution 
to  (16)  and  (17)  for  updating  the  velocity  field. 

For  a  time-invariant  pressure  gradient,  the  steady-state  solu¬ 
tions  of  (16)  and  (17)  are  given  by 


ttss 


^Giys 
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Let  Ui  and  V;  be  the  initial  velocity  components  at  time  f  , .  Then, 
by  assuming  the  pressure  gradient  remains  constant  during  the 
time  step.  At,  the  velocity  field  at  t,  +  At  is  given  by 


(u,v)\u+At  =  e 


:{Ui,Vi)  + 


-At  S 


(  'U'CC  «  Vc . 


(22) 


where  the  above  equation  is  the  analytic  solution  to  (16)  and 
(17)  evaluated  at  /.,  +  At.  We  apply  this  update  to  all  interior 
grid  nodes  to  obtain  the  velocity  field  inside  the  droplet  for  the 
current  time  step. 

The  last  piece  needed  for  updating  the  velocity  is  to  extend  it 
from  inside  the  droplet  to  outside.  In  Section  IV-C1,  the  velocity 
field  is  needed  to  update  the  level  set  function.  But  in  order  to 
do  this  properly,  it  must  be  extended  into  the  boundary  and  ex¬ 
terior  nodes  of  the  computational  domain.  This  is  accomplished 
by  letting  the  velocity  components  diffuse  into  the  exterior  re¬ 
gion  (see  Fig.  15),  which  ensures  a  continuous  velocity  field  for 
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(a) 


(b) 


Fig.  15.  Velocity  field  extension.  Illustration  showing  before  and  after  results 
of  extending  the  velocity  field:  (a)  Without  extension  and  (b)  with  extension. 

updating  the  level  set.  We  now  give  the  details  of  this  extension 
algorithm. 

First,  fix  the  values  of  uss  and  ?;ss  at  the  interior  grid  nodes 
and  set  the  edge  node  values  of  the  computational  domain  to 
be  zero.  At  each  boundary  and  exterior  grid  node  (excluding  the 
edge  nodes),  compute  a  value  for  uss  and  vss  using  this  formula: 

I  _  ^ssli+l.j  ^ss|i— l,j  “t"  ^ss|i,j+l  T  ^Ss  \i,j— 1 

where  a  similar  formula  is  used  for  vss  and  (i,j)  are  the  row 
and  column  coordinates  for  each  grid  node.  Iterate  this  process 
a  fixed  number  of  times  for  all  boundary  and  exterior  grid  nodes. 
With  each  iteration,  the  velocity  values  diffuse  further  away 
from  the  interior  region.  For  our  simulations,  we  used  50  iter¬ 
ations  to  extend  uss  and  vss,  which  are  then  used  in  (22).  This 
guarantees  that  the  velocity  field  will  always  be  extended.  Our 
method  is  nothing  more  than  Jacobi  iterations  executing  on  all 
boundary  and  exterior  nodes  and  acting  on  velocity  values.  In 
[31]  and  [48],  the  method  they  use  for  extending  the  velocity 
field  is  based  on  a  convection-type  equation  that  propagates  ve¬ 
locity  data  from  the  interior  region  into  the  boundary  and  ex¬ 
terior  nodes.  However,  we  do  not  use  their  method  because  it 
is  computationally  more  complex.  Another  technique  for  gen¬ 
erating  velocity  fields  that  also  preserves  the  distance  function 
property  of  the  level  set  function  is  given  in  [47],  but  this  is  a 


TABLE  II 

Pinch  Time  Versus  Grid  Resolution 


Grid  Resolution 

Pinch  Time  (ms) 

"o  Dev. 

84x84 

127.9 

-0.9 

90x90 

126.4 

-2.0 

96x96 

125.5 

-2.7 

102x102 

126.0 

-2.3 

108x108 

129.0 

+0.0 

114x114 

132.7 

+2.9 

120x120 

133.9 

+3.8 

130x130 

140.0 

+8.5 

140x140 

131.2 

+  1.7 

150x150 

139.8 

+8.4 

161x161 

125.0 

-3.1 

174x174 

134.2 

+4.0 

187x187 

125.4 

-2.8 

201x201 

129.8 

+0.6 

very  expensive  computation.  We  prefer  our  method  because  it 
is  simpler  and  gives  excellent  performance, 

D.  Final  Cleanup 

After  updating  the  level  set  function,  pressure,  and  velocity 
fields,  there  remain  two  final  tasks.  Reconditioning  the  level  set 
function  and  choosing  the  next  time  step.  Once  completed,  the 
program  loops  back  to  Section  IV.C  to  continue  the  simulation. 

1 )  Correct  Level  Set:  Despite  its  ingenuity,  the  level  set 
method  does  have  problems.  Since  we  are  using  the  fluid 
velocity  to  update  <j>,  it  is  highly  likely  that  the  level  set  will 
become  distorted  and  introduce  numerical  inaccuracies  [47]. 
This  requires  periodically  resetting  <f>  so  that  it  is  always  close 
to  being  a  distance  function.  This  is  done  by  explicitly  finding 
the  zero  level  set  of  (f>,  which  represents  the  droplet  boundary, 
and  recomputing  the  signed  distances  at  each  grid  point  in 
the  computational  domain.  We  speed  up  this  calculation  by 
using  a  coarse  sampling  of  the  boundary  for  computing  signed 
distances  of  grid  nodes  far  from  the  boundary.  For  closer  grid 
nodes,  we  use  a  finer  sampling.  The  advantage  of  keeping  it  a 
distance  function  is  that  it  increases  the  accuracy  of  computing 
spatial  derivatives  of  <fi.  In  addition,  it  ensures  |V</>|  «  1,  which 
increases  the  accuracy  of  computing  curvature  with  (21)  be¬ 
cause  the  denominator  is  close  to  unity.  Other  methods  exist  for 
maintaining  the  distance  function  character  of  the  level  set  (see 
[31],  [47],  [48]),  but  we  decided  to  use  a  more  straightforward 
approach. 

The  other  main  problem  with  the  level  set  method  is  that,  even 
if  it  is  updated  with  a  divergence  free  velocity  field,  it  does  not 
preserve  mass  [31].  In  general,  it  tends  to  lose  mass  as  the  simu¬ 
lation  progresses.  This  is  mainly  due  to  inherent  numerical  dif¬ 
fusion  in  the  discretization  of  (20).  We  alleviate  this  problem  by 
adding  an  appropriate  constant  offset  to  <f>  at  each  time  step.  This 
ensures  global  mass  conservation  because  the  constant  offset  af¬ 
fects  the  size  of  the  zero  level  set  (see  Fig.  14).  The  mass  is  mea¬ 
sured  by  computing  the  enclosed  area  of  the  zero  level  set,  which 
is  directly  proportional  to  the  mass  (by  incompressibility).  If 
there  is  more  than  one  droplet,  say  after  a  split,  then  different 
constants  are  added  to  the  regions  of  the  level  set  corresponding 
to  those  droplets.  Hence,  mass  is  conserved  individually  for  each 
droplet. 
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Fig.  16.  Droplet  splitting  experimental  results  with  simulation  overlay.  Six  sub-frames  showing  the  video  snapshots  of  the  experiment  (courtesy  of  CJ  Kim  at 
UCLA).  The  three  electrodes  shown  in  each  frame  have  activation  voltages  (from  left  to  right)  of  25,  0,  and  25  V.  Each  electrode  is  approximately  square  with  a 
side  length  of  1.4  mm.  The  dashed-line  droplet  outlines  are  from  Fig.  19  and  show  a  direct  comparison  between  experiment  and  the  simulation  including  contact 
angle  saturation  and  hysteresis,  (a)  Time:  0  ms.  (b)  Time:  33.33  ms.  (c)  Time:  66.67  ms.  (d)  Time:  100  ms.  (e)  Time:  133.33  ms.  (f)  Time:  166.67  ms. 


Unfortunately,  the  constant  offset  does  not  ensure  local 
mass  conservation,  such  as  in  the  pinching  region  of  a  split¬ 
ting  droplet,  which  can  cause  the  droplet  in  our  simulation  to 
‘hesitate’  while  splitting  and  sometimes  get  stuck.  The  two 
left  and  right  ends  would  bulge  and  pull  apart  as  usual,  but  as 
the  neck  joining  them  became  thin  it  stopped  moving.  This 
was  completely  erroneous  because  the  velocity  field  inside  the 
droplet  dictated  that  it  should  split  apart.  One  reason  for  this  is 
that  the  level  set  method  does  not  lose  mass  at  equal  rates  in 
different  regions  of  the  domain.  Hence,  a  constant  offset  cannot 
properly  correct  for  this.  In  addition  to  this,  the  dynamics  of 
droplet  pinching  is  not  resolved  very  well  because  the  grid 
resolution  is  fixed,  uniform,  and  quite  coarse  in  the  pinching 
region  (i.e.,  only  two  to  five  grid  points). 

Recently,  [53]  introduced  a  particle  level  set  method  that 
ensures  global  and  local  mass  conservation.  And  in  [54],  a 
method  for  adaptive  refinement  is  described  that  can  resolve 
fine-scale  dynamics.  However,  the  first  method  is  computation¬ 
ally  intensive  because  of  the  number  of  seed  particles  needed 
to  adequately  reconstruct  the  level  set  as  well  as  the  particle 
reseeding  routines  necessary  to  make  the  algorithm  work.  And 
the  second  method  leads  to  more  involved  data  structures  and 
coding.  Therefore,  we  opted  for  the  following  simpler,  faster 
method  for  correcting  the  splitting  problem. 

First,  we  check  for  potential  splitting  of  the  droplet  by  looking 
for  thin  necking  regions  in  the  flow.  This  is  done  by  using  infor- 


TABLE  III 

Simulation  Parameters 


Y-L 

Sat 

All 

Units 

L 

4.406 

4.406 

4.406 

mm 

U0 

1500 

500 

50 

mm/sec 

to 

2.94 

8.81 

88.12 

msec 

Po 

16.3 

16.3 

16.3 

N/ttr 

Re 

117.6 

39.2 

3.92 

non-dim. 

Ca 

0.01854 

0.006181 

0.000618 

non-dim. 

T 

7402.8 

2467.6 

246.8 

non-dim. 

0 

47539.1 

47539.1 

47539.1 

non-dim. 

Numbers  for  three  simulations  are  listed  here:  (Y-L)  uses 
the  Young-Lippmann  theory,  (Sat)  adds  in  saturation,  and 
(All)  includes  saturation  and  hysteresis.  Each  simulation 
uses  a  different  value  of  Uq  so  that  the  maximum  non- 
dimensional  velocity  is  close  to  unity.  This  also  causes  Re, 

Ca,  r,  and  fo  to  differ  as  well. 

mation  contained  in  the  level  set  function,  <f>,  itself.  If  it  is  not 
close  to  splitting,  then  we  do  nothing.  Otherwise,  we  modify 
4>  by  slightly  decreasing  its  height  in  a  small  region  around 
the  pinch  point  at  each  time  step.  This  prevents  the  level  set 
from  getting  stuck  and  allows  it  to  complete  pinch-off  without 
drastic  modification  to  the  level  set  function.  In  Table  II,  we 
present  simulation  results  for  the  grid  resolution  versus  time  to 
pinch-off  for  the  splitting  case  discussed  in  Section  V-A.  As  can 
be  seen  the  splitting  time  of  the  simulated  droplet  only  varies  by 
a  few  percent  from  the  108  X  108  grid  resolution  case  used  in 
Section  V. 
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(c) 


Fig.  17.  Young-Lippmann  model.  Three  simulation  frames  showing  splitting 
behavior  under  the  ideal  Young-Lippmann  theory:  (a)  time:  0.98  ms,  (b)  time: 
5.32  ms,  and  (c)  time:  10.73  ms. 

2)  Choose  Time  Step:  The  next  time  step  is  chosen  adaptively 
by  the  following  formula  [31] 

_  h 
max(|n|  +  |v|) 

where  h  is  the  grid  spacing  of  the  computational  domain,  u  and 
v  are  the  velocities  at  the  current  time  step,  and  the  maximum 
is  taken  over  all  grid  nodes  in  a  thin  band  around  the  zero  level 
set.  It  is  not  necessary  to  consider  the  whole  domain  because  we 
only  update  level  set  values  inside  the  thin  band.  This  formula  is 
based  on  the  Courant-Friedrichs-Lewy  (CFL)  condition,  which 
specifies  the  largest  time  step  that  can  be  taken  and  still  allow 
the  simulation  to  remain  stable  [55].  It  is  common  to  have  some 
safety  margin  by  choosing  a  smaller  time  step,  but  this  is  un¬ 
necessary  because  the  HJ-WENO  method  in  Section  IV-C1  has 
built-in  artificial  dissipation  which  enhances  stability.  For  more 
details,  see  [31]  and  [55]. 

V.  Results:  Theory  Versus  Experiment 

In  the  following  sections,  the  physics  of  the  splitting  droplet 
experiment  are  described  in  detail.  Next,  modeling  and  simu¬ 


(c) 


Fig.  18.  Saturation  model.  Three  frames  showing  simulation  results  when  sat¬ 
uration  is  included:  (a)  time:  2.72  ms,  (b)  time:  10.68  ms,  and  (c)  time:  17.47 


lation  results  of  the  experiment  are  presented.  We  discuss  the 
various  physical  phenomena  affecting  the  motion  of  the  droplet 
(i.e.,  saturation  and  hysteresis)  and  how  the  simulations  were 
modified  to  capture  these  effects.  Finally,  we  compare  our  sim¬ 
ulation  with  a  different  EWOD  experiment  to  show  that  a  model 
developed  for  one  type  of  experiment  is  predictive  for  a  new  ex¬ 
periment. 

A.  Theory,  Simulation,  and  Experiment  for  a  Splitting  Droplet 

In  Fig.  16,  an  overhead  view  of  an  EWOD  device  with  three 
electrodes  running  left  to  right  is  depicted  with  a  splitting 
droplet.  The  voltage  actuation,  from  left  to  right,  is  25,  0,  and 
25  V  and  is  constant  throughout  the  split.  In  frame  (a),  an  initial 
near  circular  droplet  is  shown  just  before  voltage  activation. 
After  the  voltage  is  turned  on,  the  liquid-gas  interface  over  the 
left  and  right  electrodes  deforms  and  induces  a  low  pressure 
region  there.  The  regions  where  no  voltage  is  activated  remain 
at  high  pressure.  In  the  subsequent  frames,  the  droplet  is  pulled 
from  the  left  and  right  sides,  while  it  is  pushed  in  from  the  top 
and  bottom.  The  droplet  elongates  along  the  horizontal  dimen¬ 
sion  and  is  being  pinched  in  the  vertical  direction.  This  causes 
two  smaller  droplets  to  form  on  the  left  and  right  sides,  with  a 
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Fig.  19.  Simulation  results.  Six  subframes  showing  simulation  snapshots  when  saturation  and  hysteresis  are  included  (each  frame  corresponds  to  the  experimental 
frames  in  Fig.  16):  (a)  time:  0.00  ms,  (b)  time:  33.33  ms,  (c)  time:  66.67  ms,  (d)  time:  100.00  ms,  (e)  time:  133.33  ms,  and  (f)  time:  166.67  ms. 


thin  neck  joining  them.  The  neck  eventually  gets  so  thin  that  it 
snaps  due  to  capillary  instability.  The  two  smaller  droplets  then 
continue  moving  to  the  left  and  right  electrodes  because  of  the 
pressure  differential  created  from  the  voltage  actuation.  Finally, 
the  two  droplets  come  to  rest  on  the  two  25  V  electrodes.  The 
total  time  of  this  experiment  is  approximately  167  ms. 

Next,  we  present  three  simulations  showing  the  effects  of  the 
various  physics  at  the  boundary.  The  first  simulation  is  for  con¬ 
tact  angle  variations  obeying  the  ideal  Young-Lippmann  curve 
(see  the  Y-L  curve  in  Fig.  9).  Next,  we  simulate  droplet  split¬ 
ting  assuming  only  contact  angle  saturation  (see  the  saturation 
curve  for  two  plates  in  Fig.  9).  Finally,  we  show  a  simulation 
that  includes  saturation  and  hysteresis.  In  Table  III,  we  list  the 
pertinent  parameters  of  each  simulation.  The  velocity  scale,  Uq, 
is  chosen  so  that  the  maximum  nondimensional  velocity  during 
the  simulation  is  close  to  unity. 

1 )  Ideal  Young-Lippmann:  In  Fig.  17,  we  have  a  simulation 
of  droplet  motion  when  no  contact  angle  saturation  or  hysteresis 
is  being  modeled.  As  can  be  seen,  the  general  shape  of  the  split¬ 
ting  droplet  is  not  the  same  as  in  the  experiment.  In  frame  (a), 
just  after  the  voltage  is  turned  on,  the  droplet  shape  has  much 
more  of  a  bulge  in  the  center  than  shown  in  the  experiment.  This 
becomes  more  pronounced  in  frame  (b),  with  two  thin  necks  de¬ 


veloping  between  the  three  bulging  parts  of  the  droplet.  Finally, 
in  frame  (c),  the  droplet  has  split  into  three  pieces  instead  of  two 
as  in  the  experiment.  Also,  the  total  time  to  complete  the  split 
and  reach  equilibrium  is  10.73  ms,  which  is  15.6  times  faster 
than  the  experiment. 

Since  no  saturation  or  hysteresis  is  being  modeled,  the  simu¬ 
lated  EWOD  force  is  much  larger  than  it  is  in  reality.  This  causes 
the  droplet  to  be  pulled  apart  so  fast,  that  the  middle  region  is 
never  able  to  become  a  thin  neck.  As  a  result,  three  satellite 
droplets  are  created  instead  of  two.  In  fact,  the  z  curvature  of 
the  liquid-gas  interface  (i.e.,  the  EWOD  force)  is  so  large  that 
the  x-y  curvature  component  is  practically  negligible.  This  is 
why  the  droplet  does  not  resist  being  pinched  in  two  places;  the 
EWOD  force  here  dominates  the  large  curvature  forces  induced 
by  the  pinched  regions. 

2)  Saturation:  For  the  simulation  shown  in  Fig.  18,  we  have 
added  the  effect  of  contact  angle  saturation.  The  splitting  motion 
of  the  droplet  now  looks  much  closer  to  the  experiment.  As  the 
droplet  is  pulled  apart,  a  single  thin  neck  joins  the  two  bulging 
ends.  The  neck  then  breaks,  allowing  the  two  droplets  to  come 
to  rest  on  the  left  and  right  electrodes.  However,  the  time  scale 
is  still  not  correct.  The  time  to  reach  equilibrium  here  is  17.47 
ms,  which  is  9.6  times  faster  than  the  experiment. 
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(c)  (d) 


Fig.  20.  Bulk  droplet  motion  experimental  results  with  simulation  overlay.  Four  subframes  show  video  snapshots  of  the  experiment  (courtesy  of  C.  J.  Kim  at 
UCLA).  A  time-varying  sequence  of  voltages  is  applied  to  the  eight  electrode  pattern  so  as  to  make  the  droplet  move  right,  up,  and  then  left.  Each  electrode  is 
square  with  a  side  length  of  1.4  mm.  All  device  parameters  here  are  the  same  as  for  the  splitting  experiment  shown  in  Fig.  16  except  the  electrode  pattern  is 
different.  The  dashed-line  droplet  outlines  show  a  direct  comparison  between  the  experiment  and  a  simulation  including  contact  angle  saturation  and  hysteresis 
=  0.09).  (a)  Time:  0.0  ms.  (b)  Time:  31.6  ms.  (c)  Time:  79.9  ms.  (d)  Time:  136.6  ms. 


Including  saturation  does  slow  the  droplet  down,  but  it  is  still 
not  enough.  In  reality,  the  line  pinning  induced  contact  angle 
hysteresis  is  extremely  significant  in  all  wetting  phenomena. 
Hysteresis  slows  down  motion  by  reducing  the  pressure  gradient 
(see  Section  III-B4).  Therefore,  it  is  not  surprising  that  this  ef¬ 
fect  must  be  included  to  accurately  simulate  droplet  speed.  The 
following  section  presents  our  final  simulation  of  splitting  with 
hysteresis  included. 

3)  Contact  Angle  Hysteresis:  In  Fig.  19,  we  show  six  snap¬ 
shots  of  our  simulated  splitting  droplet,  which  are  also  shown  as 
dashed-line  overlays  in  Fig.  16.  The  hysteresis  constant,  i£hys> 
is  0.09.  This  simulation  is  similar  to  the  one  in  Fig.  18,  except 
that  the  time  scale  is  now  correct.  The  simulated  droplet  now 
splits  in  the  same  amount  of  time  as  the  experiment,  as  shown 
in  Fig.  16. 

The  value  of  the  hysteresis  constant,  Khys  =  0.09,  was 
chosen  to  make  the  simulation  time  scale  match  the  experi¬ 
ment.  By  using  (12)  and  experimental  data  from  Fig.  9,  we 
estimate  the  contact  angle  deflection  due  to  hysteresis  to  be 


Ahys  =  6.4°.  In  [56],  they  give  a  value  of  about  20°  for 
sessile  drops  of  water  sliding  on  top  of  a  Teflon  surface.  This 
discrepancy  is  reasonable  since  the  droplet  size  and  geometry 
in  the  EWOD  device  is  different  than  in  [56]. 

Our  hysteresis  constant  is  also  analogous  to  the  contact  line 
friction  coefficient  in  [28],  where  they  treat  contact  line  fric¬ 
tion  as  an  extra  forcing  term  that  is  proportional  to  the  velocity 
of  the  contact  line.  In  their  case,  the  forcing  term  has  units  of 
force  per  unit  contact  line  length.  By  scaling  their  friction  force 
by  the  ratio  of  contact  line  length  to  volume  for  a  droplet  in  an 
EWOD  device  (to  put  it  into  units  of  force  per  unit  volume),  we 
can  include  this  as  a  body  force  term  in  the  Navier-Stokes  equa¬ 
tions.  After  going  through  the  same  derivation  in  Section  III-A2 
we  obtain  equations  similar  to  (6)  and  (7),  except  the  coeffi¬ 
cient  of  the  velocity  term  has  an  extra  positive  term  added  to 
1 2(  L/ H Hence,  the  coefficient  is  larger  than  before.  If  we 
ignore  the  velocity  time  derivative  term  in  (6)  and  (7),  the  extra 
friction  force  is  equivalent  to  multiplying  the  pressure  gradient 
by  a  constant  smaller  than  one  (i.e.  iThys)-  In  fact,  one  can  show 
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that  our  choice  of  hysteresis  constant  corresponds  to  a  contact 
line  friction  coefficient  of  0.5664  Newton-seconds  per  square 
meter,  which  is  comparable  to  the  values  listed  in  [28]  for  a 
column  of  fluid  comprised  of  deionized  water  between  two  pary- 
lene  coated  electrodes. 

However,  one  must  keep  in  mind  that  these  estimates  are 
rough  approximations.  The  hysteresis  constant  is  not  an  exact 
model  nor  does  it  capture  line  pinning.  But  our  goal  was  not  to 
model  line  pinning  or  contact  angle  hysteresis  in  great  detail. 
The  point  is  that  a  simple  scaling  constant  is  all  that  is  needed 
to  produce  simulations  that  approximately  capture  the  shape  and 
speed  of  droplet  motion. 

B.  Simulation  Versus  Moving  Droplet  Experiment 

In  this  section,  we  compare  our  simulation  to  another  exper¬ 
imental  case  to  give  more  supporting  evidence  for  our  model. 
The  EWOD  device  shown  here  has  eight  electrodes  arranged 
in  a  square-like  pattern.  A  predetermined  voltage  sequence  was 
used  to  actuate  the  droplet  so  that  it  moves  to  the  right  first,  then 
up,  and  finally  to  the  left.  All  constants  in  our  simulation  are  the 
same  as  those  used  in  Section  V-A3  (i.e.,  ifhys  =  0.09)  and  the 
same  computational  grid  resolution  (108  x  108)  is  used. 

From  Fig.  20,  it  is  evident  that  the  simulation  follows  the 
experiment  fairly  well.  The  match  is  not  exact,  however,  and 
this  is  mainly  because  line  pinning  is  not  taken  into  account  in 
our  model.  But  the  overall  motion  and  time  scale  is  correct. 

VI.  Conclusion 

In  this  paper,  we  have  presented  a  model  and  numerical  simu¬ 
lation  of  droplet  motion  inside  an  electrowetting  device.  Starting 
from  the  full  Navier-Stokes  equations  we  obtained  a  reduced 
order  model,  similar  to  Hele-Shaw  type  flow,  that  captures  the 
bulk  dynamic  behavior  of  EWOD  driven  droplets  in  two  dimen¬ 
sions.  The  key  part  of  our  analysis  is  including  contact  angle 
saturation  and  a  simple,  computationally  efficient  model  of  hys¬ 
teresis  in  order  to  match  the  experimental  data.  Our  simula¬ 
tion  results  show  how  these  two  physical  phenomena  affect  the 
motion  of  the  droplet.  When  all  effects  are  included,  our  sim¬ 
ulations  compare  favorably  with  the  experiments,  but  are  not 
an  exact  match.  The  main  reason  for  this  is  our  model  does 
not  include  contact  line  pinning,  which  is  observable  in  the  ex¬ 
periments.  Our  numerical  implementation  is  fast,  simple,  and 
readily  lends  itself  to  control  algorithm  design.  By  using  the 
level  set  method,  our  simulation  is  able  to  easily  capture  droplet 
splitting.  The  computing  times  of  all  simulations  (in  MATLAB) 
varied  between  three  and  six  minutes,  which  shows  the  speed  of 
our  method. 
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